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ABSTRACT: Recently, Kazakov, Gromov and Vieira applied the discrete Hirota dynamics 
to study the finite size spectra of integrable two dimensional quantum field theories. The 
method has been tested from large values of the size L down to moderate values using 
the SU(2) x SU(2) principal chiral model as a theoretical laboratory. We continue the 
numerical analysis of the proposed non-linear integral equations showing that the deep 
ultraviolet region L — > is numerically accessible. To this aim, we introduce a relaxed 
iterative algorithm for the numerical computation of the low-lying part of the spectrum 
in the U(l) sector. We discuss in details the systematic errors involved in the computation. 
When a comparison is possible, full agreement is found with previous TBA computations. 



Contents 



|l[ Introduction |] 

Y-system equations § 

p] A second formula for the ground state energy [] 

Numerical implementation | 

f3T| Discretization ^ 

372| Fixed point search 



3.3| Algorithm for the excited state 9o 
3~4| Algorithm for the excited state 6qo 



3.5| Alternative algorithm for the ground state 



||[ Numerical Results || 

PI Ground State go 

y| First excited state G 

3| Excited state 6 00 Q 

4] Summary tables and L — > limit [ll] 



IsL Conclusions Il2 



1. Introduction 

Finite size corrections to the complete spectrum of quantum field theories are an impor- 
tant issue with many applications. In particular, they are relevant to the aim of testing 
Maldacena's AdS/ CFT duality relating type IIB superstring propagating in AdS^ x S 5 
and M = 4 SYM y]. In the 't Hooft planar limit, the gauge theory composite operators 
have anomalous dimensions related to the spectrum of a finite size integrable super spin 
chain The so-called wrapping corrections are crucial to check the correspondence and 
have been the subject of major investigations in recent years [^p. 

As shown by Liischer [||], the leading corrections to the mass gap in rather general 
relativistic models can be computed in terms of their infinite size scattering matrix. In 
integrable models, the 5-matrix is one of the central objects and is known exactly. This 
suggests that the complete size corrections can be in reach. Indeed, for continuum models 
admitting an integrable discretization, it is possible to determine the spectrum of all states 
by solving the associated (nested) Bethe Ansatz equations. Their continuum limit leads 
to the Destri-de Vega non-linear integral equations (DdV NLIE) [g]. When an integrable 
discretization is not available, it is possible to determine the finite size correction to the 
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ground state energy by means of the thermodynamical Bethe Ansatz (TBA) g, §]. The 
general idea is that, in relativistic (Euclidean) 2d models, the free energy is related to 
the ground state energy in finite volume by a modular transformation exchanging spatial 
extension and inverse temperature. Usually, the method leads to an infinite system of 
coupled non linear integral equations which can be rewritten in the functional Y-system 
form which is universal, i.e. symmetry based pp. 

Apparently, the very nature of the TBA construction forbids the computation of ex- 
cited energies. Nevertheless, various extensions have been proposed [|6p. Indeed, there 
are cases under full control, like for instance the Sine-Gordon model where the very 
same Y-system form of the TBA describes excited states. They are associated to solutions 
of the Y-system with different boundary conditions and analytical properties. Due to this 
important remark, the Y-system is the most promising object for the study of the full spec- 
trum. 

In the AdS / CFT case, duality with string theory plays a crucial role by mapping the 
discrete spin chain describing finite length operators to states of a continuum cr-model. In 
the end, it allows to apply the general idea of TBA and to obtain the remarkable formu- 
lations described in [ fiO| ] 1 . The resulting Y-system leads to a very involved set of coupled 
non-linear integral equations where numerical methods are welcome, if not mandatory. 

Here, we begin a systematic investigation of such numerical algorithms by going back 
to the initial analysis of Kazakov, Gromov and Vieira who choose the Principal Chiral 



Model (PCM) as a theoretical laboratory. In [11]], the authors (GKV) assume that the Y- 
system provides a complete description of the full spectrum provided suitable solutions 
are considered. GKV describe a new systematic way to deal with the infinite component 
Y-system. They stress the important fact that the Y-system if closely related to the Hirota 
bilinear equation For finite rank symmetry group of the integrable model, the Hirota 
dynamics can be solved in terms of a finite number of functions of the spectral param- 
eter. Thus, the Y-system is successfully reduced to a finite system of non-linear integral 
equations (NLIE) rather suitable for numerical methods. 

In the GKV analysis, the proposed method is partially tested on the SU{2) x SU{2) 
Principal Chiral Model (PCM), equivalent to the 0(4) cr-model. Measuring all dimen- 
sional quantities in terms of the mass gap, the only parameter which describes the spec- 
trum is the continuous system size L. The low-lying part of the spectrum (including the 
ground state energy) is computed for 10 _1 < L < 2 by a numerical iterative algorithm. 
The general strategy starts from an approximate solution to the NLIE at large L. Then, 
the numerical solution is found by iteration of the NLIE and the size L is progressively 
reduced down to the desired value. The ground state and the first excited state energies 
can be compared with [^] where they are computed by means of coupled DdV-like equa- 
tions, apparently unrelated to the GKV NLIE 2 . The agreement is good and reasonably 



1 Actually, the string is quantized in light-cone gauge which means that relativistic invariance is broken. 
This requires to compute the free energy of a mirror model which is not the same as the original theory. 
2 Although, a priori, it is not excluded that the two formulations can be related by some change of variable. 
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compatible with the expected behaviour in the deep ultraviolet region L <C 1 which how- 
ever remains unexplored at the smallest considered size, L = 0.1. In a subsequent paper, 
the analysis of [0] has been extended to the SU(N) x SU(N) PCM [0]. For N = 3, the 
spectrum of various excitations is extracted as a function of the size L down to L = 1CT 3 . 
Unfortunately, for N = 4 already L ~ 1 is not accessible due to instabilities in the iterative 
solution of the NLIE. 

In this paper, we reconsider the GKV integral equations for the N = 2 model going 
down to very small size values providing a definitive test of their validity, at least in this 
region. We accomplish this task by a numerical analysis and present some refinements 
of the numerical strategy which, in principle, could be useful for the study of N > 2 
models. In particular we adopt a straight discretization of the NLIE which allows for 
a very clean control of the systematic errors. Also, we improve the iteration procedure 
by introducing various relaxation parameters turning out to be crucial in order to reach 
L = 1CP 8 . These very small values of the size are even smaller than those achieved in 
previous studies of the deep UV region where a comparison with asymptotically free field 
theory is possible ^ [| [l3p. The outcome of our analysis is a very good agreement between 
the GKV equations and other known results. Besides, our numerical implementation is 
very simple and easily under control. 

The plan of the paper is the following. In Sec. @, we briefly summarize the GKV 
approach and its associated NLIEs. In Sec. (||), we describe our numerical algorithms for 
the ground state and two excited states of interest. Finally, in Sec. (§|), we discuss the 
numerical results with particular attention to the role of systematic errors and relaxation 
parameters. 



2. Y-system equations 

In this section we briefly summarize the NLIE that we are going to solve numerically. 
They are derived in full details in [[H]]. Here we just sketch the main ideas and introduce 
the necessary quantities. 

The SU{2) x SU{2) Principal Chiral Model, equivalent to the 0(4) a-model, is de- 
scribed by the action 

S = J d 2 xTr(h- l d a hf , h G SU(2). (2.1) 

This model is asymptotically free. The infinite size spectrum contains a single physical 

-2-7T 

2 

particle with a dynamically generated mass m = hyy e e ° , A-uv being the ultraviolet cut- 
off. The associated state transforms in the bifundamental of SU(2) x SU{2). The elastic 
scattering matrix is known and can be written in terms of the scalar factor JT^ ] 



S {9) = 
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In general, one can consider states with N particles and polarization defined by the num- 
ber of left and right spins down 3 . In this paper we shall consider the so-called U(l) sector 
where all spins are up (ferromagnetic polarization). The infinite size energy of such states 
can be computed from the formula 

N 

E = m ^cosh^fli), (2.3) 

i=X 



where the physical rapidities {0i}i=i solve the asymptotic Bethe Ansatz [ |14| ] 

N 

e -zL sinh(^) = _ T-| 5 2(g. _ g .y ^ 

In the TBA approach, the finite size ground state energy is computed starting from the 
infinite size thermodynamics associated to the above equations. If we denote by L the size 
in units of the mass gap, then the result is the following simple formula for the energy 

E {L) = -\ ! d0cosh(7T0) log(l + Yo), (2.5) 

where Yq is obtained from the solution of the TBA equations. These are integral equations 
which, upon reasonable physical assumptions, can be written in the celebrated Y-system 
form 4 

Y+Y n = (1 + + Y n ^), neZ. (2.6) 

The physical boundary conditions predict the large 6 behaviour of Y n . All Y n with 
tend to a constant, while Yo is exponentially suppressed according to 

Yo(0) ~ e -i (2.7) 

As we explained in the introduction, the central problem is that of solving the Y-system 
since the main assumption is that it describes all states and not only the ground state. 

The GKV method to study Eq. ([D|) is based on the relation to the Hirota equation that 
can be written as (<$ is the conjugate of <3?) 

T n (e + j) T n T n m T n+l {8) = * (* + y) * (* - y) ■ ( 2 -*) 

where the functions T n are related to Y n by the relation 

YM = uy, = «pMjL$ _ i. ( , 9) 

* (9 + f J * [9 - f J * [9 + f ) * [6 - f 
Eqs. ( ^ , ^) are gauge invariant according to 

Y n {6) — > y n (fl), 



3 Here, left and right refer to the two SU(2) factors, or wings. 

4 The notation is / ± "' ± (i) = f(x ± when we have p plus or minus signs. 
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T n {6) 



in\ /„ in\ _ 



(2.10) 



The Hirota equation is integrable (it admits a discrete Lax pair) and this allows for a com- 



plete analysis of its solutions [jTTy. 



The main trick of GKV is a clever continuation from the case of large L where a fun- 
damental important simplification occurs. Indeed, the L — > oo suppression of Yq leads to 
a decoupling of the y -system in the two sets {V n } with n > or n < 0. Both can be writ- 
ten, according to ( |2.9[ ), in terms of a suitable solution of the Hirota equation and thus we 
identify two independent solutions of it. This picture can be continued to finite L as care- 
fully shown in [|ll]]. The resulting two independent solutions are no more decoupled and 
give nontrivial T n for all n on both wings of the SU(2) x SU(2) Dynkin diagram. These 
two solutions describe the same state and must be related by a gauge tranformation g{9). 
This consistency condition gives a single NLIE for the function g which is precisely the 
equation we are going to study. Once this equation is solved, the function Y" associated 
to the desired excited states is available, and the formula for the energy simply becomes 



N 



E(L) = m ^cosh(Tr^) - - / d9 cosh(7r 9) log(l + Y ). 



(2.11) 



Let us introduce the definitions 



A{6) = (g + (6)) 2 , K (9) 



1 d 
2irid9 



log Sl{9). 



(2.12) 



The function A{9) is determined by the NLIE [11] (a star denotes as usual convolution 
with respect to the rapidity argument) 



N 



A = -e- L cosh ^ e ) n Si [9 - 9i + - exp [K * log 



A- 1 

\A\ - 1 



K+ + * log 



A—l 



log 



A- 1 



where the physical rapidities are self consistently determined by 



N 



_ e zL sinh(ir 6i) tftfi - 8j) exp 



2i Im I K * log 



A - 1 



|^| - 1 ° \A\ - 1 
(2.13) 



(2.14) 



e=e % 



According to Eq. ( [2.11| ), the energy is given by 



N 



E(L) = ^cosh(7T0 t ) - I d9 cosh(7T0) log 



A-l 

\a\ - r 



(2.15) 



This formula can also be used for the ground state which is associated with N = 0. In this 
case, all factors Sq must be replaced by 1, leaving of course the kernel Kq untouched. 
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2.1 A second formula for the ground state energy 

Remarkably, an alternative NLIE for the ground state energy is also presented in [pi]]. It 
is a self -consistent equation for the function f{6) which reads 5 

2 



f{0) = 2T 1 {6) 



T [8 + - + i 



exp \-L cosh(7r/9) - 2s * log |T ++ 1 2 ] , (2.17) 



where T n are obtained from / according to 

T n -! = n + 

and where the function s(8) is 

s(9) = 



n 



7t(40 2 + n 2 ) 
1 



(2.18) 



(2.19) 



2 COSh(7T0) ' 

Once Eq. ( |2.17 ) is solved, the energy is given by Q2.5D with the following expressions for 

Tx{9)f{B) 



T (6 + h + i0 



(2.20) 



We shall discuss the numerical accuracy and efficiency of this second approach comparing 
it with the general formula for N particle states specialized at N = 0. 



3. Numerical implementation 



3.1 Discretization 



A numerical implementation of Eqs. ( |2.13 , 2.14 ) has been presented in [ ]lT| ] in terms of a 
short Mathematica code exploiting various high level functions for numerical integration 
and interpolation. Instead, our approach will be that of presenting a simple algorithm 
which can kept easily under control in order to check the continuum limit and the bias 
due to systematic errors. To this aim we cut-off the rapidity 8 in the interval \8\ < A and 
split it introducing the discrete points 

H ~ 1 ^ A£ n , n= 1.....M. (3.1) 



6n = A 



-1 + 2 



M 



For numerical efficiency, we define the fixed vectors and matrices 



5(2) 



_ e -Lco*(*B n ) Sq (e n + ~\ , n=l,...,M, 



2 log 2 



7T. 



2nA 
M — 2 



n = 0, 



n = -M, . . . ,M, 



(3.2) 



5 Note that the singular convolution is more explicitly 



f(0 = - i <% 



*J R 4(e-£ + | +l 0)= + l 7T J R 4 (fl-£ + f)» + l 



/(O-j/W- 



(2.16) 
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~(3l P' 71 = °> 

(2] ~(2) 

2im=^ m+M+1 , 71,771 = 1,...,.^, 

(31 ~(3l 

C = 2;_ m+ M + i, n,m = l,...,M, 

They are evaluated with a (large) fixed number of digits of order 200. The issue of numer- 
ical precision will be discussed later. We now present the numerical algorithms for the 
ground state, the first excited state which has a zero Bethe root, and another state which 
is non trivial from the Bethe Ansatz point of view having two opposite Bethe roots. 

3.2 Fixed point search 

The NLIE that we want to solve takes the form of a fixed point equation f(x) = x. It 
can be solved by iteration when the map a;n>/(a;)isa contraction. In other words, one 
simply sets x n+ \ = f(x n ) and the sequence {x n } converges to the solution x* as n — > oo, 
independent on the starting point. In actual cases, one needs to start in a suitable small 
neighborhood of x* . Then, the map can be linearized in terms of the small differences 
<5n = x n — x* and one has to iterate 

Sn+i = a6 n , a = f'(x*). (3.3) 

Convergence is exponentially good for \a\ < 1, otherwise the sequence diverges. Follow- 
ing Jacobi and Seidel, a simple general trick to deal with the case \a\ > 1 is to introduce a 
relaxation parameter and consider some average between the old and new values of the <5 
variables. In other words, one considers the following modified equation 

x n+1 = (1 - X)x n + Xf(x n ), (3.4) 

where A is a real number. The advantage is clear. Around the fixed point one has 

S n+1 = {l- X + Xa)6 n . (3.5) 

This means that in the unstable case, \a\ > 1, we can enforce stability by simply choosing 

2 

a > 1 : < A < 0, 

a ~ 1 2 (3-6) 

a < -1 : < A < 



1-a' 



In the p-component case, where x n G C p , the same argument based on the linearization 
around the fixed point, suggest to replace A by a p x p matrix and to apply (^) to each 
eigenvalue of the matrix V/. More involved schemes can be devised replacing ( |3.5[ ) by 
higher order relaxation processes as in the inertial relaxation discussed in |TJ] or also 



finding adaptive algorithms for the step by step choice of A Q17|]. In this paper, we shall 



discuss the simplest implementation of the relaxation idea as illustrated above. 
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3.3 Algorithm for the excited state 9q 



Let us begin our discussion with the 1-particle state 60 . It is discussed in |Tl|] and is the 
lightest excited state. So, it enters the mass gap computation. By symmetry, there is only 
one Bethe roots which is fixed at zero. This simplifies a lot the computation. 

(k) 

The function A(6) becomes a sequence of vectors {A n }k>o,i<n<M which takes the 
initial value 

A® = Z n l) = -e-^Orfn) 5o (g n + (37) 
Then, the main loop starts. Defining 

r« = log^^-, (3.8) 

we have 

A^ = (1 - \) A k - 1] + (3.9) 



+A Z n 1] exp 



2A N 



M -2 



/ w m \ n nm ' m ^nm' m j ' n 



The parameter A is a relaxation parameter introduced according to the previous discus- 
sion. It is crucial at small L as we shall discuss. The algorithm in [ |TT| ] is obtained for A = 1 
and works for L of order 10 _1 . We shall push the computation several orders of magni- 
tude beyond this value. In this deep UV region, the above map is no more a contraction 
and convergence is lost. Nevertheless, a moderate value of A £ (0, 1) 6 will turn out to be 
enough to restore convergence. The weights w m are used to approximate integrations by 
summations. They can be taken at the rough value w m = 1 or improved as in Simpson or 
higher order integrations. For our purposes this is not a crucial issue as we shall explain 
later. 

The very same algorithm can be applied to the ground state. The only necessary 
change is the replacement of all So factors by unity. In conclusion the parameters entering 
the numerical algorithm are 

A, M, A. (3.10) 



3.4 Algorithm for the excited state 600 

This is a 2-particle state. By symmetry, it is associated with a pair of opposite rapidities 
±6* solving the Bethe Ansatz equations. Again, we start from the method presented in 
[pl|], but add the important relaxation parameters. Indeed, it is convenient to choose 
two independent relaxation times for the function A and for the Bethe root(s). This is 
a fundamental point otherwise convergence is soon lost. 

6 We never observed the necessity of over-relaxation. 
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Let us denote by {9 a } a =i,2/ with B 2 = — #i , the two Bethe rapidities. They satisfy 



the BA equations (actually it is one equation due to the parity symmetry !) pi]], 
L^t^-ijy^S^ rt 1] ) =0. 

6=1 

(3.11) 

The initial value of the Bethe roots is not critical and a possible choice is just to solve the 
above equation suppressing the last term in the left hand side. The modified main loop is 
then 

a£ ] = (1 - M) + (3.12) 



x exp 



x 

N 



— V w (Z {2) r ik ~ 1] - Z {Z) r {k ~ 1] \ -r (k - 1] 



M - 2 

m=l 



In addition, we must provide an update rule for the Bethe roots. We propose 

t k) = (1 - A 2 ) 0a k ~ 1] + A 2 da k ~ 1] . (3.13) 

In summary the parameters entering the numerical algorithm are 

A, M, Ax, A 2 . (3.14) 

As a remark, we note that for more complicated states with more Bethe roots, we still have 
two relaxation parameters, one for the function A and one for the Bethe roots. 

3.5 Alternative algorithm for the ground state 

The numerical implementation of the solution of ( 2.17| ) is completely analogous to that 



described in Sec. (f3J|). The initial value of the discretized / function is 
4. Numerical Results 

We present here our numerical results aimed at understanding (a) the systematic errors 
associated with finite values of A, M, (b) the role of the relaxation parameters A^. Let 
us begin with the cut-off A. In principle, one should consider a fixed A and send M — > 
oo in order to reach the continuum limit. The result is a function of A which has to be 
extrapolated at A — > oo. In practice, as soon as A is large enough to cover the support of the 
function A(9), one observes independence on A with exponentially increasing accuracy. 
In the following we shall set 

n A , /81og 10\ 

A = arccosh f — | — J . (4.1) 
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When tia = 1/ w e have the exponential factor e~ L cosh ( 7r A ) = io~ 8 which is a reasonable 
definition of a small quantity in this problem. We shall work with several hundreds of 
digits which will be enough in all cases where convergence is achieved. 



4.1 Ground State 



We compared the computation of the ground state energy based on fl2.13| ) or fl2.17p . The 
first equation provides much more accurate results and shall be discussed in this section. 
In the final summary table we shall present results for the second method too. As a pre- 



liminary step, we start with L = lCP 1 which is the smallest value considered in [11]. We 
present the time history of the energy iterates at relaxation A = \ in Fig. (^). The dis- 
cretization is M = 350. One sees that there is easy convergence to a value which has a 
mild residual dependence on A. This dependence can be observed together with the de- 
pendence on M in Fig. @. A natural simple guess is to assume that for large enough A 
(all considered cases should be all right from this point of view), the relevant parameter is 
the density of points A/M. This hypothesis is tested in Fig. @. One sees that indeed the 
dependence on A and M is under control. The considered values of A are already asymp- 
totic and the dependence on M is quite accurately linear 7 . As an interesting additional 
information, we provide in Fig. @ the real and imaginary parts of the function A{6). One 
sees the initial profile before the iteration loop as well as the final equilibrium value. 

The same analysis can be repeated at the very smaller value L = 1CT 6 . Again, we 
present the time history of the energy iterates at relaxation A = A in Fig. (||). The dis- 
cretization is M = 500. The dependence on A and M is shown in Fig. @. The scaling 
plot showing dependence on A/M is Fig. (0). Finally, we provide in Fig. (§) the real and 
imaginary parts of the function A{6). 

This small value of L permits to emphasize the role of relaxation. This is shown 
in Fig. @ where we show the time history of the energy as the relaxation parameter is 
reduced from 1, where instability is observed, to 1/3 where convergence is all right. 



4.2 First excited state 9o 

We present our numerical data for the excited state 6o following the same scheme as for 
the ground state. In particular, we show for L = 10 _1 the time history of the energy 
iterates at relaxation A = \ in Fig. (10). The discretization is M = 350. The dependence on 
A and M is shown in Fig. (|TT1 ). The scaling plot showing dependence on A/M is Fig. (|T^). 
Similar plots for L = 10~ 6 can be found in Figs. (|l3|, |l4[ |l5|)). Finally, we provide in Fig. dlrj) 
the real and imaginary parts of the function A{6). 

Again, at L = 10~ 6 , we can show the role of relaxation. This is shown in Fig. ( |T7[ ) 
where we show the time history of the energy as the relaxation parameter is reduced from 
1, where instability is observed, to 1/2 where convergence is all right. 



7 In principle one can improve the convolutions by improved numerical integrations. Actually, this would 
be an incomplete improvement and we have checked that residual linear corrections 1 /M do remain. 
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4.3 Excited state 6 



oo 



This is the next-to-lowest excited state, at least for small enough size L. The relaxation 
algorithm depends now on two independent parameters: Ai for the update of the function 
A and A2 for the update of the Bethe root(s). We present in Figs. ( |T8| , |19| ) the time history of 
the energy of the state 9oo as well as that of 6\ for M = 50, = 2, and fixed Ai = 1/2. As 
the second parameter A2 is reduced, we move from an oscillating regime to a convergent 
one. The onset of oscillations instead of an exponential instability is compatible with the 
fact that we are dealing with a two component system. 

Next, using the parameters shown in Tab. we present at L = 1CT 6 similar plots 
to those that we have discussed for the other states. In particular, with discretization 
M = 600, we show the dependence on A and M in Fig. (|20|), the scaling plot in Fig. 
and in Fig. p2[ ) the real and imaginary parts of the function A{9). 



Finally, in Fig. (23) we show the profile of A{9) which is obtained after convergence 
at various sizes L. These could help, at least in principle, in the formulation of suitable 
proposals for the analytic size dependence of this important function. 

4.4 Summary tables and L — > limit 

Our final results, obtained after extrapolating to M — > 00, are summarized in three ta- 
bles. For the ground state, they are shown in Tab. (|l|). The first column reports the results 
obtained with the TBA NLIE of [0] down to L = 10~ 6 . The next column shows the re- 
sults obtained using the equation d2.17| ). Finally, the third column shows the results which 
we obtain using Eq. ( 2.13| ) which turns out to be much more efficient and accurate. The 



Y-system results have been obtained reducing the size by two orders of magnitude com- 
pared with [^]. The similar Tab. (||) presents our results for the energy of the state 60 while 
we present in Tab. @ the energy of 9 00 for which there are no available results obtained 
with other methods. 

One can check that there is a very good agreement showing that the GKV equations 
are working perfectly in the very small size limit. A better precision could be achieved 
by simply increasing M in order to reduce the effect of the subleading correction. To the 
aim of testing the GKV equations, we honestly believe that the quality of our result is 
convincing. 

In Fig. (p4|), we show the plot of the three energies as functions of L. The mass gap is 
clearly reproduced given the agreement between Tab. ([l]) and the analysis of [^]. There is 
an additional check that we can perform on our data, and that follows the analysis of [O]. 
Indeed, all the three considered energies should have the limit E — > —1/4 as L — > 0. 
This limit is definitely out of reach for the numerics of [|Tl|]. Fitting our numbers in the 
range 10 ~ 8 < L < 10~ 5 by means of a (naive) quadratic polynomial in 1/ log L we obtain 
the three estimates 

ground state : -0.249(3), 9 : -0.2477(1), 6 00 : -0.2439(3). (4.2) 

These extrapolations are rather close to the predicted limit —0.25. The quoted errors are 
simply those inherited from the data. There is no attempt to estimate the systematic error 
due to subleading corrections in L — > which are apparently larger for the excited states. 
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L 


E (NLIE) [0] 


#o (Y-System) Eq. ((2.17) 


E (Y-System) 


Parameters 


10 1 


-11.273364587(1) 


-11.273(1) 


-11.2734(2) 


M = 200 - 350, A = 1/2 


10 2 


-127.22634373(1) 


-127.22(8) 


-127.226(2) 


M = 250 - 400, A = 1/2 


1 n — 3 
10 


— 1343.4090793(1) 


— 1343(3) 


10/10 A C ( 0\ 

— 1343.46(2) 


M = 300 — 450, A = 1/2 


io- 4 


-13865.238816(1) 


-1.386(4) x 10 4 


-13865.8(2) 


M = 350 - 500, A = 1/2 


IO" 5 


-141563.8217(1) 


-1.415(5) x 10 5 


-141564(5) 


M = 350 - 500, A = 1/10 


IO" 6 


-1436683.423(1) 


-1.434(7) x 10 6 


-1.4367(1) x 10 6 


M = 350 - 500, A = 1/10 


io- 7 






-1.4526(2) x 10 7 


M = 350 - 500, A = 1/15 


io- 8 






-1.4651(4) x 10 8 


M = 350 - 500, A = 1/15 



Table 1: Ground states energies. The first column is taken from [|^]. The second column reports our 
best results. The last column shows the range of discretization M and the relaxation parameter. 



L 


E(& ) (NLIE) [0] 


E{& ) (Y-System) 


Parameters 


IO" 1 


-3.004108884(1) 


-3.004109(8) 


M = 200 - 350, A = 1/2 


IO" 2 


-69.83802786(1) 


-69.8380(3) 


M = 350 - 400, A = 1/2 


10~ 3 


-901.2815867(1) 


-901.282(9) 


M = 400 -450, A = 1/2 


io- 4 


-10260.214298(1) 


-10260.2(1) 


M = 450 - 500, A = 1/2 


IO" 5 


-111091.0324(1) 


-111091(2) 


M = 500 - 550, A = 1/2 


IO" 6 


-1172575.496(1) 


-1.17258(2) x 10 6 


M = 550 - 600, A = 1/2 


io- 7 




-1.21947(1) x 10 7 


M = 550 - 600, A = 1/2 


io- 8 




-1.25637(1) x 10 8 


M = 550 - 600, A = 1/2 



Table 2: 0o energies. The first column is taken from [|^]. The second column reports our best 
results. The last column shows the range of discretization M and the relaxation parameter. 



L 


E{e 00 ) (Y-System) 


Parameters 


IO" 1 


12.423491(8) 


M = 200 - 350, Ai = 1/2, A 2 = Ai/5 


IO" 2 


52.4704(3) 


M = 200 - 350, Ai = 1/2, A 2 = Ai/5 


10~ 3 


123.324(3) 


M = 250 - 400, Ai = 1/2, A 2 = Ai/5 


io- 4 


-1515.85(2) 


M = 250 - 400, Ai = 1/2, A 2 = Ai/5 


IO" 5 


-35284(1) 


M = 300 - 450, Ai = 1/2, A 2 = Ai/5 


IO" 6 


-5.0557(1) x 10 5 


M = 300 - 450, Ai = 1/5, A 2 = Ai/10 


io- 7 


-6.2486(2) x 10 6 


M = 200 - 450, Ai = 1/5, A 2 = Ai/20 


io- 8 


-7.2033(6) x 10 7 


M = 200 - 350, Ai = 1/5, A 2 = Ai/20 



Table 3: 0oo energies. The first column reports our best results. The second column shows the 
range of discretization M and the relaxation parameters. 

5. Conclusions 

The recent proposal of GKV [ |TT| ] provides a quite general method to compute finite size 
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correction to the full spectrum of two dimensional integrable models. It involves non lin- 
ear integral equations that can be treated in the full space of physical parameters only by 
means of numerical methods. Whenever integrable discretizations are not available, the 
calculation of excited levels is based on certain assumptions. For this reason, it is impor- 
tant to compare different methods as well as achieve accurate numerical predictions. In 
this paper we have worked out the small size limit of the SU(2) x SU(2) Principal Chiral 
Model. To this aim, we have tested a numerical implementation of the NLIE of [|TT|]. We 
have explored the possibility of solving them by iteration in the case of the ground state 
and of two additional excited states. We found that small L values require to introduce 
relaxation constants in order to achieve convergence. This has to be done independently 
for the NLIE and Bethe root evolution. We hope that this investigation will be useful in 
the analysis of the SU{N) x SU(N) Principal Chiral Model for general N along the lines 
of [jT2|] . We believe that this first steps are necessary in order to attack the full Y-system 
equations for the AdS / CFT problem [[jl]] having all the systematic errors of the numerics 
under control. 
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Figure 1: L = 10 _1 , M = 350. Time history of the ground state energy iterates with relaxation 
A = 1/2. The notation of the legend is A = n A x • • • . The three values of n A are used to check 
cutoff independence. 
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Figure 2: L = 10 1 . Ground state energy after convergence, as a function of 1/M, the discretiza- 
tion roughness, for three values of the cutoff A. The notation of the legend is A = tia x • • • . 
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Figure 3: L = 10 1 . Ground state energy, after convergence, as a function of A/M, the discretiza- 
tion density, for three values of the cutoff A. The notation of the legend is A = n A x • • • . 
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Figure 4: L = 10 1 . Profile of A{6) for the ground state, before and after convergence. 
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Figure 5: L = 10~ 6 , M = 500. Time history of the ground state energy iterates with relaxation 
A = 1/10. The notation of the legend is A = n A x • • • . The three values of n A are used to check 
cutoff independence. 



-1.405x10 6 



-1.410x10 6 



-1.415x10 6 



-1.420x10 6 



-•- A = 2x 

- A - 3 x 
-♦- A = 4 x 



-1.425X10 6 



• 1 

0.0020 0.0022 0.0024 0.0026 0.0028 0.0030/w 



Figure 6: L = 10 6 . Ground state energy after convergence, as a function of 1/M, the discretiza- 
tion roughness, for three values of the cutoff A. The notation of the legend is A = n A x • • • . 
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Figure 7: L = 10 6 . Ground state energy, after convergence, as a function of A/M, the discretiza- 
tion density, for three values of the cutoff A. The notation of the legend is A = n A x • • • . 
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Figure 8: L = 10 6 . Profile of A(9) for the ground state, before and after convergence. 
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Figure 9: Convergence of the energy of the Qq state at L = 10~ 6 , M = 100, and ua = 2, for various 
values of the relaxation parameter A. The computation is done with 250 digits. When the plots 
start oscillating wildly, convergence is lost and the numerical accuracy rapidly decreases. The 
plateau which is observed in the first phase of the evolution at A = 7/10 is in agreement with the 
convergence at A = 1/2. 
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Figure 10: L = 10 _1 , M = 350. Time history of the 8 state energy iterates with relaxation 
A = 1/2. The notation of the legend is A = n A x • • • . The three values of n A are used to check 
cutoff independence. 
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Figure 11: L = 10 1 . 6 state energy after convergence, as a function of 1/M, the discretization 
roughness, for three values of the cutoff A. The notation of the legend is A = tia x ■ ■ ■ . 



-20- 



E(©o) 
-300800 - * 

-301 000 - 

-301 200 - 



0.008 0.010 0.012 0.014 0.016 0.018 0.020 M 



Figure 12: L = 10 1 . O state energy, after convergence, as a function of A/M, the discretization 
density, for three values of the cutoff A. The notation of the legend is A = tia x ■ ■ ■ . 
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Figure 13: L = 10~ 6 , M = 600. Time history of the 8 state energy iterates with relaxation 
A = 1/2. The notation of the legend is A = tia x ■ ■ ■ . The three values of tia are used to check 
cutoff independence. 
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Figure 14: L = 10 6 . 60 state energy, after convergence, as a function of 1/M, the discretization 
roughness, for three values of the cutoff A. The notation of the legend is A = tia x ■ ■ ■ . 
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Figure 15: L = 10 6 . 8 state energy, after convergence, as a function of A/M, the discretization 
density, for three values of the cutoff A. The notation of the legend is A = tia x ■ ■ ■ . 
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Figure 16: L = 10 6 . Profile of j4(0) for the Q state, before and after convergence. 
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Figure 17: Convergence of the energy of the 8 state at L = 1CT 6 , M = 100, and ua = 2, for 
various values of the relaxation parameter A. The computation is done with 250 digits. When the 
plots start oscillating wildly, convergence is lost and the numerical accuracy rapidly decreases. 
The plateau which is observed in the first phase of the evolution at A = 7/10 is in agreement with 
the convergence at A = 1/2. 



-23- 



11 11 




Figure 18: Convergence of the energy of the G o state at L = 10 6 , M = 50, and n A = 2, for Ai = | 
and various values of the relaxation parameter A 2 . 
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Figure 19: Convergence of the Bethe root &\ at L = 10 6 , M = 50, and n A = 2, for Ai = | and 
various values of the relaxation parameter A 2 . 
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Figure 20: L = 10 6 . G o state energy, after convergence, as a function of 1/M, the discretization 
roughness, for three values of the cutoff A. The notation of the legend is A = tia x ■ ■ ■ . 
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Figure 21: L = 10 6 . 8 00 state energy, after convergence, as a function of A/M, the discretization 
density, for three values of the cutoff A. The notation of the legend is A = tia x ■ ■ ■ . 
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Figure 22: L = 10 6 . Profile of A{8) for the G o state, before and after convergence. 
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Figure 23: Profile of A{9) for the 9 o state after convergence at various sizes L. 
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Figure 24: Summary plot showing the size dependence of the energies of the three considered 
states. 



